An algorithm of calculating transport parameters of thermoelectric materials using single Kane band model with Riemann integral methods

We develop an algorithm called SKBcal to conveniently calculate within minutes the thermoelectric transport parameters such as reduced Fermi level (η), electronic thermal conductivity (κe), lattice thermal conductivity (κl), Hall factor (A), effective mass (m*), quality factor (β) and theoretical zT within the framework of single Kane band (SKB) model. The generalized Fermi–Dirac integral for SKB model is integrated by left Riemann integral method. A concept of significant digits of relative error is involved to determine the accuracy of calculation. Furthermore, a combined program of "For" and "While" is coded to set up an iteration for refining the reduced Fermi level. To easily obtain the quality factor, we re-derive the expression into a formula related to carrier mobility. The results calculated by SKBcal are consistent with the data reported in the literatures.

One can analyze transport properties of a given thermoelectric material using the single parabolic band (SPB) for conduction in a single parabolic band 1 or the single Kane band (SKB) model for conduction in a single nonparabolic band. Either model can be used to calculate parameters of the reduced Fermi level (η), electronic thermal conductivity (κ e ), density-of-states effective mass at the band edge (m 0 * ) and other transport parameters by using experimentally determined thermopower and Hall concentration. One could therefore use the carrier concentration (n) dependence of the quality factor (β), which is calculated using the above models as a strategy to optimize the dimensionless figure of merit zT = σS 2 /κ, where σ, S, and κ are the electrical conductivity, thermopower, and thermal conductivity, respectively.
For the energy band theory in solids, the crystal momentum dependence of the energy would be nonparabolic if the direct energy gap is small at the appropriate band edge or the carrier concentration is high. Kane found a highly nonparabolic band in InSb and good agreement between experimental data and the calculated results on fundamental absorption and its dependence on n-type impurity concentration 2 .
To begin with, the reduced Fermi energy is first calculated by using experimentally determined thermopower (Seebeck coefficient) and two Fermi-Dirac integrals. The Lorenz number is then calculated by three Fermi-Dirac integrals using the as-obtained reduced Fermi energy. As seen in the following equations 3,4 ,

OPEN
Department of Physics, National Changhua University of Education, Changhua 500, Taiwan. * email: liucj@cc.ncue. edu.tw The Hall mobility (μ H ) is given by Eq. (12), where is the reduced Planck constant, e the elementary charge, C l a combination of some elastic moduli relevant to the acoustical vibrations 3,8 , m b * the average band mass (density-of-states effective mass for each band pocket or valley) 4,6 , m I * the inertial mass of the carriers along the conducting direction (inertial effective mass) 4,9 , and E def the deformation potential coefficient 3,4,[6][7][8] . Moreover, the relation between density-of-states effective mass for each band pocket and m 0 * is given as 4 , where the N v is the number of degenerate carrier pockets 4,6,7 . Therefore, Eq. (12) can be rewritten as the following, Finally, the theoretical zT can be estimated using Eq. (15), x n x + εx 2 m (1 + 2εx) 2 + 2 k/2 dx, where β is the quality factor of thermoelectrics and can be expressed as Eq. (16), Chasmar and Stratton first described that the maximum zT could be determined by the quality factor 10 . As seen in Eq. (16), it can be readily seen that the combination of the lager number of band valley (N v ), small inertial effective mass (m I * ), and low lattice thermal conductivity (κ l ) is beneficial to obtain high-performance thermoelectric materials. It can also be used as one of the strategies to select or engineer a good thermoelectric material.
The SKB model has been adopted to estimate the material parameters and optimize the thermoelectric performance for several material systems such as InSb, PbTe, PbSe, PbS, SnTe, SnS, Cu 2 SnSe 4 , ZrNiSn, and elemental Te [2][3][4]6,7,[11][12][13][14][15][16][17][18] , which indicates the importance of the SKB model in understanding thermoelectric materials with nonparabolic model, especially for chalcogenide materials. However, the knowledge of band gap (E g ), the band degeneracy (N v ), and the ratio of longitudinal to transverse effective mass (K) of the material are required to calculate the thermoelectric parameters using the SKB model.
The Fermi-Dirac integral for the SKB model seems a little complex than that of the SPB model and difficult to compute; it cannot be easily carried out by some mathematical software. However, the calculation can be realized using approximate values with mathematical methods. In this work, we calculate the Fermi-Dirac integral with a simple mathematical method of Riemann integral using software programs written in Python. Furthermore, we calculate the parameters of the reduced Fermi level (η), Lorenz number (L), electronic thermal conductivity (κ e ), lattice thermal conductivity (κ l ), Hall factor (A), carrier concentration (n), carrier mobility (μ), densityof-states effective mass (m 0 * ) at the band edge, quality factor (β), and dimensionless figure of merit (zT), which are then saved into a file with the CSV file format. In Riemann method, the accuracy can be determined by the partition size for each subinterval. The relation between the partition size and the significant digits of relative error (SDORE) of some calculated parameters would be discussed.

Mathematics
Derivation of generalized Fermi-Dirac integral for SKB model. Zawadzki  where n L m k is a general formula of GFDI, f the Fermi distribution, ε = k B T/E g the reciprocal reduced band gap, and x = E/k B T the reduced energy. In the model of calculating transport properties, the relaxation time must be introduced in the integral function. Thus, Eq. (17) should be rewritten as a complete GFDI for the SKB model given by 3 , where n R m k is the complete GFDI for the SKB model, f the Fermi distribution, ε = k B T/E g the reciprocal reduced band gap, and x = E/k B T the reduced energy. Besides, τ (x) is the relaxation time for the acoustic phonon scattering and can be expressed as where D (x) is the density of state as a function of x. In general, the density of states in SKB model is considered as a function of energy and is given by 3 , where m 0 * is the density-of-states effective mass at the band edge, the reduced Plank constant, E g the band gap, and E the energy. In order to combine both the variables of relaxation time and density of state, the definitions of the reciprocal reduced band gap (ε = k B T/E g ) and the reduced energy (x = E/k B T) should be substituted into the expression of the density of state. Thus, the density of states can be rewritten as www.nature.com/scientificreports/ Now, the relaxation time can be further rewritten as After simplification, we would obtain a simplified formula given by In fact, the total relaxation time consists of two components: acoustic phonon scattering and optical phonon scattering. The combined formula is expressed as In the calculation of transport parameters using the SKBcal algorithm, we assume the acoustic phonon scattering is dominant within the framework of the SKB model. Therefore, the symbol of τ (x) would represent the relaxation time for scattering by acoustic phonons or the total relaxation time in this work. Now, we can substitute the expression of relaxation time into Eq. (18) and obtain There is a constant term preceding the integral and is only dependent on the index k. Thus, we can ignore the constant term in some parameter calculations, such as the reduced Fermi level, Lorenz number, Hall factor, and carrier concentration. However, the expression of Hall mobility has two GFDIs of one with the index k = -2 in numerator and the other with k = 0 in denominator. For convenience, the constant term in Hall mobility was taken out of the numerator part in GFDI, which is already shown in Eq. (14). Then, the form of GFDI with the SKB model in Eq. (26) can be adopted to calculate all of the parameters in this study. Here we repeat the Eq. (26) for a convenient reading, Derivation of transport parameters within the framework of SKB model. In this section, our derivation of the formulas for transport parameters will involve a simplified Kane's band model 20 , which is also called first order nonparabolicity 21 . As seen in the shape function, the nonparabolicity of energy band is determined by the index k. The band structure is parabolic if k is 0. In order to describe the non-parabolic band model, only the first two terms in Eq. (27) are used, where γ 1 is unity due to the parabolic band. Substituting Eq. (28) into the definition of density-of-states effective mass, where γ 2 is the reciprocal band gap. One can also obtain this relation by using a simplified Kane's model given by, Equation (30) shows the dependence of the wave number on the energy. Besides, the effective mass (m*) entering into the transport properties of carries for a spherical energy band of arbitrary shape has the energy dependence of where ε = k B T/E g is the reciprocal reduced band gap and x = E/k B T the reduced energy. Furthermore, the expression of the carrier mobility is given by where τ (x) is the relaxation time for the acoustic scattering process. The average of carrier mobility is then expressed as As compared to Eq. (14), it can be seen that the inertial effective mass (m I * ) is equivalent to the average band mass (m b * ). The contributions of scattering by charged impurity centers and polar optical phonons are ignored in the calculation within the framework of the SKB model. Besides, we do not consider the anisotropic materials. Hence, the equality between the two effective masses might arise from only considering the acoustic scattering events and isotropic material. Moreover, the Hall factor, A, which is used to convert the Hall mobility and Hall concentration into charge carrier mobility, is the product of anisotropy factor (A k ) and statistical factor (A τ ) and is given by 3 Finally, the dimensionless figure of merit zT is expressed as With the above derivation for n, μ, S, and L, we can rewrite zT as www.nature.com/scientificreports/ where B is a constant and is the quality factor. According to the position of B in the formula, one can readily see that large B leads to high zT values. The quality factor reported by literatures is usually given by 4,6,14 which is one third of B. Despite of that, the quality factor is a scale to estimate the performance of thermoelectric materials. In our SKBcal coding, we adopt β as the quality factor for consistency with literatures.

Code implementation Generalized Fermi-Dirac integral (GFDI) within the framework of SKB model. According to the expressions of transport parameters, the GFDIs of
0 should be calculated. However, it is a difficult task to obtain these values by direct integration. Hence, we decide to solve the GFDIs by a simple numerical method of Riemann integral. In the beginning, the differentiation of the Fermi distribution is derived as By substituting Eq. (44) into Eq. (3), the GFDI can be expressed as where η is the reduced Fermi level. It is well known that the size of regular partition (dx) would affect the accuracy of the resulting value in the Riemann method. In order to get the idea of their correlation, we calculate the significant digits of relative error (SDORE) for each GFDI using different partition sizes. The relative error is obtained via dividing absolute error by GFDI integration value. The absolute error is usually estimated using the mean value theorem and is given by where M is the greatest first derivative over the interval [a, b], p the magnitude of regular partition, and f (x) represents the function of each GFDI. However, the way of predicting the maximum error is not quite suitable to GFDIs within the framework of the SKB model. The slope on the curve of f (x) is sometimes very large and even becomes infinite. Therefore, we adopt a simple method to estimate the error. As seen in Fig. 1 The predicted error using the above method is better than that estimated by the mean value theorem because of the slope problem of GFDIs in the SKB model.
In order to get the idea of the influence of regular partition on the error, we calculate the significant digits of relative error (SDORE) of each GFDI using different partition sizes. The SDORE in this study is given by www.nature.com/scientificreports/ which shows the level of accuracy in calculations. As seen in Fig. 2, the SDORE increases significantly with smaller regular partition. For simplicity, the band gap and absolute temperature are set as 0.05 eV and 300 K, respectively. The SDORE of each GFDI are almost at a similar level in the same regular partition. Overall, the SDORE shows little dependence of on the reduced Fermi level. However, the SDORE for 0 F 1 −2 and 0 F 1/2 −4 does show some increases when the reduced Fermi level turns from negative to positive. It is worth mentioning that the SDORE is about 4 if the regular partition is 0.0001, this accuracy is sufficient to perform the calculation. Hence, we would select this partition size of 0.0001 for calculating all the transport parameters.
For a complete analysis, we calculate the SDORE using different magnitudes of band gap and absolute temperature. Figure 3 shows the reduced Fermi level dependence of SDORE for different values of band gap ranging from 0.025 to 1 eV. Even though the SDORE differs in various band gaps, it is relatively small and approaches zero at high reduced Fermi level.

Reduced Fermi level (η).
In this section, we would discuss the relation between the thermopower and reduced Fermi level. The interval of the reduced Fermi level should be determined in the calculation. The reduced Fermi level can be refined and determined by experimentally measured thermopower using Eq. (1), which is rewritten as In the refinement, the method of iteration is adopted to obtain the refined reduced Fermi level since the reduced Fermi level involves with GFDIs. In order to confirm the suitable upper limit in GFDIs, the reduced Fermi level dependence of thermopower is calculated and plotted, which is shown in Fig. 4. We take the regular partition of 0.0001, the band gap of 0.05 eV, and the temperature of 300 K for the calculation. It can be readily seen that all the reduced Fermi level is less than 12 for thermopower in the range of 1 to 1000 μV/K. It is noted that materials with a small thermopower less than 10 μV/K can be viewed as metals, which are not suitable for calculation using the SKB model. However, the interval of [0, 80] is chosen in our calculation using SKBcal. Furthermore, the reduced Fermi level is about 56.1288 under the calculation conditions: thermopower = 0, partition size = 0.0001, band gap = 0.05 eV, and T = 300 K. However, we will not consider the case of thermopower = 0 in SKBcal, in particular for thermoelectric materials.
In the coding, we use the functions of "For" and "While" to build the combination of iteration. After entering the experimentally determined value of thermopower in SKBcal, a test value of zero is assigned to the initial reduced Fermi level, which is used to plugged into GFDIs. The values of GFDI integration are easily calculated by "For" loop in SKBcal. A rough reduced Fermi level is then obtained by Eq. (50). The reduced Fermi level is further refined by the While" loop. The reduced Fermi level is determined and recorded and saved in the file once the absolute value of the difference between the rough and calculated reduced Fermi level smaller than 0.0001, namely, "abs(RFL-testRFL) in the "While" loop meets the criteria of "abs(RFL-testRFL) < 0.0001". As such we can refine the reduced Fermi level to the fourth decimal place. It should be noted in our SKBcal coding that the regular partition size is set to 0.01 for the initial guess of the rough value of reduced Fermi level to speed up the runtime of iteration. After the rough value is obtained, the regular partition size is set to 0.0001 for the subsequent refinement of the reduced Fermi level and subsequent calculation of transport parameters.
(49) SDORE = − log 10 E total GFDI integration value , www.nature.com/scientificreports/ Lorenz number (L), electronic thermal conductivity (κ e ), and lattice thermal conductivity (κ l ). After the refined reduced Fermi level is obtained, the Lorenz number is then calculated immediately by Eq. (2). Likewise, the electronic thermal conductivity is estimated by the Wiedemann-Franz law using Eq. (4). The lattice thermal conductivity is then computed by subtracting the electronic thermal conductivity from the measured total thermal conductivity.
Hall factor (A), carrier concentration (n), carrier mobility (μ), and effective mass (m * ). Utilizing the SKB model, the electronic transport parameters including the Hall factor, carrier concentration, and densityof-states effective mass at the band edge can be calculated by using Eqs. (8), (10), and (11), respectively. Besides, the carrier mobility is given by The ratio of the density-of-states effective mass at the band edge relative to the rest electron mass is directly calculated in the SKBcal.

Quality factor (β) and calculated thermoelectric figure of merit (zT). Calculating and predicting
the zT is one of the important functions in our SKBcal, which is beneficial to optimizing the thermoelectric properties of a material. However, the quality factor β is one of the parameters in the theoretical zT calculation. www.nature.com/scientificreports/ As seen in Eq. (16), β involves the elastic moduli (C l ), deformation potential coefficient (E def ), and the inertial mass (m I * ). In general, there is difficulty in finding these three material parameters from literatures. However, we find an easy way to obtain β from the carrier mobility. By comparing Eqs. (12), (16), and (50), the expression of β can be rewritten as where the density-of-states effective mass at the band edge and the lattice thermal conductivity can be obtained from the previous section. As such we can calculate β and zT without knowing the values of C l , E def , and m I * in Eq. (16). www.nature.com/scientificreports/

Code availability
The algorithmfor calculating the transport properties of thermoelectric materials can be accessed through the Supplementary Information.

Data verification
To confirm the accuracy of our SKBcal coding, we extract electrical resistivity (ρ), thermopower (S), and Hall concentration (n H ) from graph data in Ref. 11 using WebPlotDigitizer 23 . We then calculate the density-of-states effective mass using our algorithm of SKBcal by setting the value of the band gap, band degeneracy, and K being 0.43 eV, 1, and 1, respectively. Table 1 shows the comparison between our calculated results (shaded area in the table) of density-of-states effective mass and the lattice thermal conductivity (κ l ) and those calculated based on the SKB model in Ref. 11 . The percent error shown in parenthesis is given by The percent error in Table 1 is lower than 4%, which should be acceptable because transport parameters used in calculation show temperature discrepancy due to different experiments (resistivity, thermopower, and Hall effect). Besides, there are two cells entered as "None" because the total thermal conductivities at 425 and 525 K are not recorded in Ref. 11 .
To further verify the SKBcal, we repeat the same procedure as above and make a comparison between calculated data obtained using SKBcal with one more set of data reported in Ref. 24 , and show the results in Table 2.   www.nature.com/scientificreports/ The value of the band gap, band degeneracy, and K is set to 0.39 eV, 1, and 1, respectively. Below the temperature of 600 K, the percent error is lower or near 5%. However, the percent error of the lattice thermal conductivity suddenly increases to 15% at 700 K. This large percent error can be ascribed to the phase transition of CdSnAs 2 at high temperatures 25 . According to both of the slopes in dρ/dT and dS/dT change around 600 K reported in Ref. 24 . This indicates that the electronic structure of Cd 0.95 Zn 0.05 SnAs 2 around that temperature has turned to a different state from room temperature. Therefore, the band gap of 0.39 eV adopted in calculation below 600 K is no longer good for T ≥ 600 K. Table 3 shows the comparison between calculated results using SKBcal and transport parameters reported in Ref. 26 . The transport parameters of Lorenz number (L), electronic thermal conductivity (κ e ), and lattice thermal conductivity (κ l ) in Ref. 26 are calculated by Mathematica using Riemann integral method. As compared to the results obtained using our SKBcal by Python, the difference between them is almost negligible.

Conclusion
We develop an algorithm called SKBcal coded by Python 3.7 to calculate and predict the transport parameters for materials with a single nonparabolic band structure within the framework of the SKB model. The generalized Fermi-Dirac integrals (GFDIs) in this study are calculated by the left Riemann integral method. To quickly obtain the reduced Fermi level, the iteration process is implemented by two stages: the first stage to obtain a rough value, which is then used at the 2nd stage for determining the refined reduced Fermi level. The determined reduced Fermi level is then used for subsequent calculation of all transport parameters. Accuracy of numerical integration is determined by the significant digits of relative error (SDORE) which is larger than 4 when setting the regular partition smaller than 0.0001. As compared to the data reported in the literature, the results calculated using SKBcal are in good agreement with the data reported by literature.